Fluid Vesicles with Viscous Membranes in Shear Flow 
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The effect of membrane viscosity on the dynamics of vesicles in shear flow is studied. We present 
a new simulation technique, which combines three-dimensional multi-particle collision dynamics for 
the solvent with a dynamically-triangulated membrane model. Vesicles are found to transit from 
steady tank-treading to unsteady tumbling motion with increasing membrane viscosity. Depending 
on the reduced volume and membrane viscosity, shear can induce both discocyte-to-prolate and 
prolate-to-discocyte transformations. This dynamical behavior can be understood from a simplified 
model. 

PACS numbers: 87.16.Dg, 47.55.-t, 87.17.Aa 



The dynamical behavior of vesicles — closed lipid 
membranes in aqueous solution — under shear flow is 
an important subject not only of fundamental research 
but also in medical applications. For example, in mi- 
crocirculation, the deformation of red blood cells reduces 
the flow resistance of microvessels. In diseases such as 
diabetes mellitus and sickle cell anemia, red blood cells 
have reduced deformability and often block microvascu- 
lar flow. Although red blood cells do not have a nucleus 
and other intracellular organelles, they are of more com- 
plex than simple lipid vesicles, since their plasma mem- 
brane has an attached spectrin network, which modifies 
its elastic and rheological properties. The deformability 
of cells and vesicles is determined by their shape, the vis- 
cosity of the internal fluid, and the viscoelasticity of the 
membrane 

The dynamical behavior of vesicles in shear flow has 
been studied experimentally 0, Q, theoretically 0, 0, 
and numerically d, 0, 0- The vesicle shape is deter- 
mined by the competition of the curvature elasticity of 
the membrane, the constraints of constant volume V and 
constant surface area S, and the external hydrodynamic 
forces. One of the difficulties in theoretical studies of 
the hydrodynamic effects on the vesicle dynamics is the 
boundary condition for the embedding fluid on the vesi- 
cle surface, which changes its shape dynamically. In some 
previous studies, a fluid vesicle was therefore modeled as 
an ellipsoid with fixed shape Q • More recently, the time- 
evolution of the shape was studied numerically using a 
boundary integral method in three spatial dimensions Q 
or an advected-field method in two spatial dimensions 
. The red blood cell membrane has also been modeled 
as an elastic capsule of discocyte shape Q. 

Two types of dynamics have been found in these stud- 
ies, a steady state with a tank-treading motion of the 
membrane and a finite inclination angle with the flow 
direction, and an unsteady state with a tumbling (flip- 
ping) motion. A transition from tank-treading to tum- 
bling motion with an increasing viscosity of the inter- 
nal fluid has been predicted for fluid vesicles with fixed 
ellipsoidal shape in three dimensions Q, and with the 
advected-field method in two dimensions 0. When the 



shape is relaxed dynamically in three dimensions, all dis- 
cocyte vesicles were surprisingly found to transform into 
prolates in shear flow, even for the smallest shear rates 
studied [5j . 

In this letter, we focus on the effect of the membrane 
viscosity on the dynamics of vesicles in shear flow. This 
is an important question, because the membrane of red 
blood cells, for example, becomes more viscous on aging 
0, 0] or in diabetes mellitus Q . Experiments indicate 
that the energy dissipation in the membrane is larger 
than that inside a red blood cell 0. Furthermore, it has 
been shown recently that vesicles can not only be made 
from l ipid bilayers, but also from bilayers of block copoly- 
mers njj. These "polymersomes" have been shown to 
have a membrane viscosity which is several orders of mag- 
nitude larger than for liposomes [ll| . 

Several mesoscopic simulation techniques for fluid flow 
have been developed in recent years. We present here the 
first simulation studies for a combination of mesoscopic 
model for the solvent and a coarse-grained, dynamically- 
triangulated surface model for the membrane. This ap- 
proach has four main advantages: (i) The membrane is 
described explicitly, so that their properties like the vis- 
cosity can be varied easily; (ii) thermal fluctuation of 
both the solvent and the membrane are fully and consis- 
tently taken into account; (iii) the method can easily be 
generalized to more complex flow geometries; and (iv) no 
numerical instabilities can occur. 

We employ a par ticle-based hydrodynamics method 
E 13 111 III EU3 El to simulate the solvent, which is 
called multi-particle collision dynamics (MPCD) [l7Hli| 
or stochastic rotation dynamics 0, 0|. This method 
was applied, for example, to flow around a solid object 
[l7j and to polymer dynamics [l^l Il8j | . The fluids in the 
interior and exterior of the vesicle are taken to be the 
same, in particular to have the same viscosity T]q. 

As the MPCD model is described in detail in Refs. 0, 
E fl5| , we can be very brief in explaining the meso- 
scopic simulation technique. The solvent is described by 
N s point-like particles of mass m s moving in a rectangu- 
lar box of size L x x L y x L z . The algorithm consists of 
alternating streaming and collision steps. In the stream- 
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ing step, the particles move ballistically and the position 
of each particle is updated according to 



T i (t + h)=T i (jt)+V i (t)h, 



(1) 



where Vj is the velocity of particle i and h is the time 
interval between collisions. In the collision step, the par- 
ticles are sorted into cubic cells of lattice constant a. The 
collision step consists of a stochastic rotation of the rel- 
ative velocities of each particle in a cell, 
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(2) 



where v cm is the velocity of the center of mass of all par- 
ticles in the cell. The matrix Q(p) rotates velocities by a 
fixed angle <p around an axis, which is chosen randomly 
for each cell. In our simulation, the angle ip — tt/2 is em- 
ployed. We apply a random-shift procedure 14] before 
each collision step to ensure Galilean invariance. 

For the membrane, we employ a dynamically- 
triangulated surface model in which the membrane is 
described by iV m b vertices which are connected by tethers 
to form a triangular network. The vertices have excluded 
volume and mass m m b. The shapes and fluctuations of 
the membrane are controlled by curvature elasticity with 
the energy |2cj . 



C 2 ) 2 dS, 



(3) 



where k is the bending modulus, and C\ and Ci are the 
principal curvatures at each point of the membrane. The 
curvature energy is discretized as described in Ref. [2l| . 
To model the fluidity of the membrane, tethers can be 
flipped between the two possible diagonals of two adja- 
cent triangles. These bond flips provide also a convenient 
way to vary the membrane viscosity ?7 m b, because it in- 
creases with decreasing bond-flip rate. In contrast to pre- 
vious studies of dynamically triangulated surfaces, which 
were all done by Monte Carlo simulations, we introduce 
a smooth bond-interaction potential, which makes the 
model amenable for molecular dynamics simulations. 

The solvent particles interact with the membrane in 
two ways. First, the membrane vertices are included in 
the MPCD collision procedure [ljj. Second, the solvent 
particles are scattered elastically or via bounce-back from 
membrane triangles. We use here the procedure sug- 
gested in Ref. |l6j for a spherical particle. 

To induce a shear flow, we employed Lee-Edwards 
boundary condition [l5l which gives a linear flow 
profile (v x ,v v ,v s ) = (7*, 0,0) in the MPCD fluid. The 
particle density was set to p = 10r7i s /a 3 (N s = 450000, 
L x = 50a, L y — L z — 30a). In experimental con- 
ditions of red blood cells and liposomes, the Reynolds 
number Re= ^pR\jr\Q is very small (Re ~ 10~ 3 ), where 
i?o = yJS/A-K is the effective vesicle radius. Therefore, 
we chose a short mean free path hyfk B T/m s = 0.025a, 
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FIG. 1: (Color online) Dependence of the average inclina- 
tion angle (0) (—tv/2 < 9 < 7r/2) on the membrane viscosity 
77m b for reduced shear rate 7* = 0.92 and various reduced 
volumes V*. The error bars are estimated from three inde- 
pendent runs. Squares and circles represent discocyte and 
prolate vesicles at V* = 0.59, respectively. Triangles and dia- 
monds represent prolate vesicles at V* = 0.78 and V = 0.91, 
where the prolate is the only stable shape. The solid lines 
and broken line are calculated by K-S theory with prolate 
(V* = 0.59, 0.78, and 0.91) and oblate ellipsoids (V* = 0.59), 
respectively. 



where k B is the Boltzmann constant and T is the tem- 
perature 0]. Then the viscosity of solvent fluid is rjo = 
20.1y/m s k B T/a 2 [3. We used k = 20k B T, N mh = 500, 
and m m b = 10m s . The volume V and surface area 
S = 405a 2 of a vesicle are kept constant to about 1% 
accuracy. With these parameters, we obtain a Reynolds 
number Re~ 0.1. The results are conveniently expressed 
in terms of dimensionlcss variables: the reduced volume 
V* = V/(4ttRq/3), the intrinsic time scale r = tioRq/k, 
the reduced shear rate 7* = 7T, and the relative mem- 
brane viscosity rj^ h = ^mb/^o^o- Details of the numeri- 
cal scheme will be published elsewhere j^] ■ 

At rj* b = 0, a vesicle exhibits tank-treading motion for 
all simulated reduced volumes in the range 0.59 < V* < 
0.97. We calculated the average inclination angles (9), 
and found them to agree very well with those obtained 
by the boundary integral method ||. 

With increasing membrane viscosity 77* b , the inclina- 
tion angle 9 decreases, as shown in Fig.^ The qualitative 
features of the simulation data are reproduced very well 
by the theory of Keller and Skalak (K-S) H0. Note that 
there are no adjustable parameters. Due to the approxi- 
mations in the K-S theory, an agreement on a quantita- 
tive level cannot be expected: (i) an ellipsoidal shape is 
assumed, which is only mimics the real shapes of vesicles, 
(ii) the flow on the surface of the droplet is not locally 
area conserving, as it must be for a incompressible mem- 
brane, and (iii) thermal fluctuations are ignored in the 
theory, but are present in the simulations. In the K-S 
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FIG. 2: Free-energy profile F(a) of the asphericity q for 
V* = 0.59 in the absence of shear flow. The sliced snapshots 
of stable (discocyte) and metastable (prolate and stomatcyte) 
shapes are also shown. 




theory, the vesicle transits from tank-treading to tum- 
bling motion when the angle 9 reaches 0. In contrast, 
we observe tumbling intermittently to occur already for 
nonzero (#), since our simulation includes thermal fluc- 
tuation. For example, the vesicle with V* = 0.78 starts 
tumbling at r]* b — 1.22. This intermittent tumbling 
smoothes out the decrease in (9) around the transition 
point, compare Fig. ^ 

We now focus on the case V* — 0.59. At this re- 
duced volume, the oblate discocyte shape is stable and 
the prolate and stomatocyte shapes are metastable in the 
absence of shear flow. Fig.|2shows the free energy F as a 
function of the asphericity a, calculated with a version of 
the generalized-ensemble Monte Carlo method 24] . The 
asphericity a = [(A 1 -A 2 ) 2 + (A 2 -A 3 ) 2 + (A 3 -A 1 ) 2 ]/[2i?4] ; 

with the eigenvalues Ai < A 2 < A3 of the moment-of- 
inertia tensor and the squared radius of gyration i? 2 = 
A1 + A2 + A3, is a convenient measure to distinguish oblate 
and prolate shapes, where a = for spheres, a = 1 for 
thin rods, and a — 0.25 for thin discs [25j. The shear flow 
changes this stability. For membrane viscosity 77* lb = 
and shear rates 7* > 1.66, the discocyte state is found to 
be destabilized and to transform into a prolate, in agree- 
ment with the results of Ref. However, for a smaller- 
shear rate of 7* = 0.92, the discocyte vesicle retains its 
shape. Speculations about shear to be a singular pertur- 
bation [5j can therefore be ruled out. 

The inclination angle 9 of prolates decreases faster than 
that of discocytes with increasing r^ b , see Fig. U At a 
large membrane viscosity of ry^ b = 1.62, the prolate en- 
ters the tumbling phase, while the discocyte remains in 
the tank-treading phase. The reason is that the disco- 
cyte has a flat dimple region and is less affected by the 
membrane viscosity than the prolate. Remarkably, for 
small shear rates 7* < 1.84, the (metastable) prolate 
starts tumbling, but after a 7r or 27r rotation, transforms 
into a tank-treading discocyte, see Fig- EI Only for larger 
shear rates of 7* > 2.76, the tumbling continues — ac- 



FIG. 3: (Color online) Time development of (a) asphericity 
a and (b) inclination angle 9, for V* = 0.59, rj^ b = 1.62 and 
7* = 1.84. The broken lines are obtained from Eqs. (J1J and 
© with Ca = 100, A = 15, and B(a) = 1.1 - 0.18a. 
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FIG. 4: (Color online) Time development of (a) aspheric- 
ity a and (b) inclination angle 9 at 7* = 2.76. The other 
parameters are set to same values in Fig. [3] 

companied by shape oscillations between prolate and dis- 
cocyte (Fig. |IJ. At intermediate membranes viscosities, 
■q* mh = 0.49 and ri* mh = 0.87, and shear rate 7* = 0.92, the 
prolate transforms into a discocyte after tank-treading 
motion for a time of (70 ± 40)r or (40 ± 20)r by ther- 
mal fluctuation, respectively. For a larger shear rate of 
7* = 1.84, the tumbling continues intermittently. 

K-S theory 0, g explains the 77^-dependence of the 
stability of tank-treading (compare Fig. QJ, but cannot 
be applied to describe the dynamics, including morpho- 
logical changes. We suggest the simple equations 

Cad = -K^dF/da + Aj* sin(26») (4) 
<9 = 0.57*{-l + 5(a)cos(26>)}. (5) 

The force OF/ da is calculated from the free energy F(a) 
of Fig. El The second term of Eq. (@J is the deformation 
force due to the shear flow. Its angular dependence can 
be derived from the shape equations of Ref. [5j , while the 
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amplitude is assumed to be independent of the aspheric- 
ity a (to leading order). Eq. (J5J is adopted from K-S 
theory 0,0. Here, B is a constant which depends on 
viscosities and ellipsoid shape. For B > 1, a steady an- 
gle 9 — 0.5 arccos(l/-B) exists and tank-treading motion 
occurs, while for B < 1, there is no stable angle and tum- 
bling motion occurs. In our case, the vesicle shape can 
be time dependent, so that B is no longer constant. For 
simplicity, we assume a linear dependence of B on the 
asphericity, B(a) — Bq — B\a. To obtain tank-treading 
discocytes and tumbling prolate, we need -6(0.2) > 1 
and S(0.7) < 1, respectively. Then, Eqs. (0} and © 
reproduce the simulated dynamics very well, see Figs. |3 
and The vesicle is found, for example, to relax after 
some tumbling to a stable, tank-treading discocyte state 
at 7* = 1.84, and to relax to a limit-cycle oscillation at 
7* = 2.76. 

The vesicle deformation due to shear depends on the 
inclination angle 9. Shear flow increases the elongation of 
a vesicle for < 9 < ir/2 (where 7* sin(26>) > 0), and the 
angles 9 of tank-treading motion belong to this region. 
On the other hand, shear flow reduces the elongation for 
— it/2 < 9 < (where 7* sin(26>) < 0) during tumbling. 
The force in the former case induces the discocyte-to- 
prolate transformation, in the latter case the prolate- 
to-discocyte transformation. With increasing membrane 
viscosity 77^, the inclination angle 9 of a tank-treading 
discocyte decreases, and larger shear rates 7* are neces- 
sary to generate the required clongational forces to in- 
duce a discocyte-to-prolate transition. 

It is also interesting to compare the effect of membrane 
viscosity ?7 m b and internal viscosity rj- m . In both cases, an 
increase of the viscosity induces a decrease of the incli- 
nation angle 9 and a transition from tank-treading to 
tumbling. However, the effect of internal viscosity rj- m is 
less dependent on the vesicle morphology. K-S theory 
H shows that the tank-treading phase of oblate vesicle 
is destabilized a little faster than that of prolate with 
an increase of r]i n (for ?7* lb = 0); the transition viscosity 
is r) in /j] — 2.8 and 3.3 of oblate and prolate ellipsoids 
at V* = 0.59, respectively. Thus, only a sufficiently 
high membrane viscosity r)* b can induce a prolate-to- 
discocyte transformation. 

The membrane viscosity of human red blood cells is 
estimated from the analysis of the tank-treading motion 
to be rjmb = 10 _7 Ns/m [4j, while a micropipette recovery- 
time technique gives ?y m b = 10~ 6 Ns/m 8]. When the 
viscosity of the external fluid is set to the same value of 
the intracellular fluid, 770 = 10 _2 Pa s, and i?o = 3.3/jm, 
the relative membrane viscosity is then found to be in the 
range rj* b = 1...10. Thus, the effect of this membrane 
viscosity is sufficiently large (compare Fig.^) to strongly 
affect the dynamics of erythrocytes. The viscoelasticity 
of membrane can be changed by varying the chemical 
composition of the solvent |lj. It is difficult to separate 
the effects of viscosity and elasticity, however. On the 



other hand, the membrane viscosity of polymersome can 
be changed by varying the polymer chain length. Thus, 
polymersomes seem to be very well suited to study the 
effect of membrane viscosity experimentally. 

To summarize, we have applied the MPCD method to a 
vesicle with viscous membrane under a simple shear flow. 
The membrane viscosity qualitatively changes the vesi- 
cle dynamics. The shear induces a discocyte-to-prolate or 
prolate-to-discocyte transformation at low or high mem- 
brane viscosity, respectively. The deformations of other 
vesicles shapes, such as stomatocyte and budded vesicles, 
will be interesting subjects in the further studies. 
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